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Abstract. We present and derive a technique for the introduction of defects 
into a beam model based on the Cosserat theory of rods. The technique is designed 
for the derivation of component models of non-ideal rods for use in MEMS devices. 
We also present a worked through example of blob/nick defects (where the rod 
has an area with an excess/lack of material) and a guide for a model with random 
pits and blobs along the length of the beam. Finally we present a component 
level model of a beam with a defect and compare it to results from a Finite 
Element Analysis simulation. We test the Cosserat model for two cases without 
any defect and four with a defect. Results are in good agreement with a maximum 
0.5% difference for the ideal case and under 1% differences for all but one of the 
defective cases, the exception being a 2% error in an extreme case for which 
the model is expected to break down. Overall, the Cosserat model with and 
without defects provides an accurate way of modelling long slender beams. In 
addition, simulation times are greatly reduced through this approach and further 
development for both component level models as well as as FEA components is 
important for practical yet accurate modelling of MEMS both for prediction and 
comparison. 
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1. Introduction 

Many MEMS devices rely on flexible rods of Si (or other materials) to act as a 'spring' 
for the device. As such, accurate models of such rods have been developed [1, 2, 3, 4] 
to facilitate simulation of these devices. Although these techniques can be quite 
sophisticated, so far work has concentrated on perfect (prismatic) rods i.e. those 
whose cross-sections do not change along the length and which have uniform material 
properties. While this is quite adequate for most common MEMS, such as the Tang 
Resonator where the system is designed to eliminate any minor defects caused by 
the manufacturing process, some more advanced devices may be quite susceptible to 
defects as may present designs if miniaturisation is to take place. 

Presently, most simulation of such devices would be carried out using Finite 
Element Analysis through a package such as ANSYS®} or FELT [5]. These 
simulations can be extremely time consuming both at the design level and during 
the calculation. Once defects are added it could be expected that the time required 
for design and to a lesser extent computation will increase dramatically. 

Another approach to modelling is component level (or network) modelling [6, 7]. 
Here each component of a MEMS such as beams, shuttles, comb drives etc. arc 
modelled individually and joined together. The resulting system requires finding the 
solution of a system of differential-algebraic equations (DAEs) [8] which is a relatively 
simple task. As such, component level modelling is very efficient for the modelling of 
new devices. Another aspect of component models is that they can also be introduced 
into FEA calculations to speed computation of certain 'regions' or 'components' of the 
system. Traditionally, beam models such as those of Euler-Bernoulli or Timoshcnko 
[9] are used in these circumstances. 

In this paper we propose a technique to model rods which builds on a recent 
model of beams based on Cosserat Theory [10] to introduce non-ideal properties such 
as defects caused by the manufacturing process. To demonstrate the viability of this 
method we provide two fully worked through and tested examples being a rod with a 
blob of excess material (or a 'nick' which is essentially the same problem) and a rod 
in which the manufacturing process has left a random jitter i.e. a series of minute 
pits and blobs with a known statistical distribution. We also run a series of tests on 
the blob/nick model where FEA simulations are used as a benchmark of the model 
developed in this paper. 

The technique is designed to be easily implcmcntablc in a component level 
simulation, allowing it to be used in rapid modelling of present and future MEMS 
devices. The design of a 'defective component' can be carried out using this approach 
and installed into a preexisting model. 

2. Preliminaries 

In the modelling of flexible rods, the Cosserat approach is a good method for both 
analytic and numerical treatments. Through a kinematic assumption that the cross- 
sections of the rod change their orientation and position but not their shape along 
the length of the rod (physically reasonable in most cases) it is possible to write a 
Lagrangian for a rod as follows [10] (henceforth we use bold-italics for matrices and 

t ANSYS®is a registered trademark of ANSYS, Inc. 
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vectors with capitals for matrices and lowercase for vectors) 



V = —(u — u)K(u — u) + (u — u)T(v — v) + -(v — v)J(v — v) ds 
Jo 2 2 

C = T- V 



T = / -Ad t r 2 + w x A\d\ ■ d t r + -wlw ds 
Jo 2 2 




(1) 



where d s r = u, d s di = u x di, dtdi = w x di, K = Kijdi ® dj, J = Jijdi ® dj, 
T = Tijdi^dj and I = Iijdi^dj. t is time and s is a label for the position of a cross- 
section along the rod length (defined to run from to L so that As — 1 corresponds 
to one unit in the rest shape), u and v define the shape of the rod when no force is 
applied anywhere (reference configuration). We use Einstein summation convention 
where Greek letters are 1 or 2 and Roman letters range from 1 to 3 for the rest of this 
paper and repeated indices indicate sums eg. :— x iUi-% 

Finding the stationary point of the Lagrangian through the usual method with 
two variables yields the following coupled equations (three dimensions each leading to 
a six dimensional, second order linear PDE) 



where / and I are external forces and torques and n(s,t) = d v V(s,t), m(s,t) = 
d u V(s, t), q = A y d 7 and h = Iw. 

As in [3] we will adopt a quasi-static assumption. This means that we assume 
that all time-dependent terms in (2) (i.e. the RHS) are set to zero but the boundary 
conditions vary with time. In the context of MEMS modelling this assumption can 
be justified be realising that the frequencies of the entire, multi-component system, 
are significantly lower than the vibrational frequency of each component. This means 
that the wavelength of eigen modes is greater than half the length of the rods. With 
this assumption we reduce (2) to a six-dimensional, second order ODE which is much 
more amenable to analytic work than the complete PDE. 

3. Material Perturbations 

In order to introduce a defect we can introduce changes to any (or all) of: the shear 
and elastic moduli G and E; the shape of the cross-section changing A, A\ and J; 
and the reference configuration through changes to u and v. Our technique involves 
introducing any of these changes as a perturbation to the ideal case. That is, we find 
a closest possible ideal case and treat differences as small. We can then expand (to 
increasingly higher orders if we like) the Cosserat equation in terms of the perturbation 
coefficient (which we denote T) of the defect. The expansion in terms of T should then 
approach the complete solution as we increase the number of terms. 

Noticing that any changes to the material will only affect n and m, we may 
rewrite (2) as (to first-order in T) n(v,u) = n^°'(v,u) + Tn^^v^u) and m(v 7 u) = 
m^(v,u) + Tm^°\v,u) and then solve the Cosserat equations as a perturbation to 
the known system with n = and m = m^. 

§ The tensor quantities above take the following values: = Sx-,G,Kx^ = K-^x = 0, K33 = E, 

J\-f = ^\ 1 EIx 1 ,J\3 = J:i\ = 0, J33 = GI33, T A7 = T33 = 0,Ti 3 = —-5^2^23 = _BAi,T 3 i = 
GVl2,T32 = — GAi and I Xl = 5a 7 ^W I\3 = h\ = 0, i"33 = Here A is the cross-sectional 

area, Ax its first mass moments and Ax 1 are its second mass moments. 



d s n + f = pAd tt r + d tt q 

d 8 m + (d s r) xn + l = qx (d tt r) + d t h 



(2) 
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Keeping only the first-order terms and maintaining the quasi-static assumption 
we can now write 

d s n {1) = 

(3) 



where 



n (i) = K^Av^ + KWAvW + [A(xW),if< ' 



m« = j(0) Atl (i) +J (i) Att (o) + \ A (xW), jH Au< ' + t( 1 W°) 



(4) 



and Av^ = - Aw« = - u« and A(x^) = #( O U(0«)#(°) T 

:= R(<p^) = ® ej). Here the commutator term with A{x^) takes into 
account the effect of the change of basis set to first order. 

This is then solved according to the boundary conditions r^(0) = r^(L) = 
and (1) (O) = <t) (1) {L) = 0. 



3.1. Series expansion in e 

As in [3] let us make a first-order expansion in e (the perturbation coefficient of the 
boundary conditions) and substitute it into (4). We can now write 

A«(°> =eAvW,AuW = e Au^ 

At/ 1 ) = Av^ + eAv^lAu^ = Au^ + eAu^ 

Here the first upper index refers to powers of T and the second to powers of e and 
v (o,o) _ ^ yy e can j us tify the exclusion of terms above 0(e) through the following 
argument: the expression we use for the ideal case is accurate to 0(e 3 ); our expression 
is accurate to 0(Te); assuming manufacturing to be quite precise then T ~ 0(e 2 ); 
thus the total accuracy for the non- ideal case ~ 0(e 3 ) 

Under these assumptions, equations (3) and (4) become d s n^ 1,0 ^ = 0, c^n^ 1 ' 1 ) = 

and 

d s m^ + A(v^)n^ + A( V ( 1 '°))n(°'°) = 
d s m^ + A(vW)n^ + A(vW)n^ 

+ A(vW)nW + A(v^)nW = 



where 



n (M) = K (o,o) Av (iS) + K (ifi) Av (o,i) + Ati (o,i) T (i,o) 

+ \a(xW),kW] AvW + ^A(x^),K^] AvW 

m (i,i) = j(o,o) Ati (i,i) + j(M) Am (o,D + T (i,o) At> (o,i) 

+ \A(xW), J(°>°)1 Au^ + \a(x^), J(°'°)l A«(°-°). 



(5) 
(6) 



4. Examples 

^J. 'BZofes' and Wicfo' 

Sometimes a rod might be created with a bump of extra material or a nick through 
errors in the manufacturing process. This would locally affect the mass moments of 
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the rod so that we would have to use A(s), A 1 (s) and A lX (s) to calculate K, J and 
T. 

Treating this defect as a small perturbation from an otherwise uniform cross- 
section we can write 

K = K^+TK^°\s), J = j(°>°) + rj( 1,0 )(s), T = TTW{s) 

where = K?.' n0) df ' 0) ® df' 0) and i^ 1 ' ) = K&' >(*)d< ' 0) ® df 0) while 

K£' 0) (s) = i^ 0) dist(s) (same for J and T). Note that T< < ) = as we choose 
the path of the ideal rod to follow the centre-of-mass of the fixed cross-sections. 

Let us find a solution to (3) where the rod deviates slightly from its reference 
frame (to order e) due to displacements and rotations of its end points. We substitute 

A„(°.°) =eu< ' 1 ),Au< ' > =euW 
& v M) =ev (i,i) tAu (i,o) =eu (i,i) 

into (5) and (6) (here the terms involving commutators vanishes as A(x^' ^) is by 
definition O(e) and must therefore be multiplied by Ai/ 0,0 ) = 0) as follows 

n (M) = #"(0,0)^(1,1) _ 0(1,1)) + ^"(1,0)^(0,1) _ 0(0,1)) + T (1,0)T U (0,1) 
m (l,l) = J(0,0) tt (l,l) + j(l,0) u (0,l) + r (l,0) ( „(0,l) _ 0(0,1)) 

which may then be used to calculate t/ 1 ' 1 ) and u^ 1 ' 1 ). 

After some work we can obtain expressions for (p^ 1 ' 1 ^ and r^ 1 ' 1 ) were we have 
made an assumption that the shape of the perturbations distribution is governed by 
the same function of s for each of K, J and T (this is completely true for nicks and 
blobs). The solutions are thus given by 

j(o,o) (i,i) = fe a,o )s _ A3fe (i,o)^ 

—J' (fc#°>l - Ask^s) - T'[K^]- 1 k^ o n 
K (o,o) r (i,i) = K (o,o) A3 J s (i,i) ds / + k (i fi ) s _ K > k (o,o)l 

-T' T [J^-\k%°ri-A 3 k^§) 

where f(s) = f(s')dist(s')ds' . We can calculate fcl 1 ' ^ and (fcn°'°^ an d fcm'°^ 

arc known from the ideal case) by ensuring that <fP~ ,1 ^{L') = and r^^iL) = 0. The 
full analytic solutions are far too complicated for inclusion here although the MAPLE 
source to generate them is available on request. 

Figure 1 shows a visual demonstration of the model were we have calculated the 
bending properties of a beam of certain dimensions and boundary conditions both with 
and without a defect. The difference is, as expected, greatest around the position of 
the defect itself. Both the change in shape, and potential energy from the defect itself 
will affect the restoring force of the beam. 

4-2. Random jitter 

The treatment of random jitter is all but identical to that of a localised defect such as 
a blob or nick. The only difference is that the distribution function dist(s) will have 
different expectation values which depend on the type and magnitude of the noise (eg. 
we may have 1 = J Q L dist(s)lds = whereas s = k). Substituting these values into the 
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Figure 1. Model of a beam demonstrating its shape when it contains/does not 
contain a defect. The blue represents the ideal case while the green beam has a 
defect located at the marked point. 



final expressions (which we have not included here) we can calculate the behaviour 
of the rod. Although beyond the scope of this paper, the same model will also be 
applicable to the calculation of noise 'distributions'. 

This second example demonstrates the power of the technique as the same 
equations can be used for two different types of physical impurities. Other impurities 
can be modelled in similar ways by changing the 'noisy' variable. 

5. Implementation 

One of the main advantages of the method outlined in [4] is that it gives us a means 
to convert the internal structure of a rod into a generalised force dependent only on 
the value of the end points (position and orientation of each end). This allows us to 
easily generate a component model of the rod for use in a component level simulator. 

Using MAPLE we can follow the same procedure to derive analytic expressions 
for the change to the effective spring matrix of the end points (a 12 x 12 matrix) 
caused by a defect. We then add our new term to the third-order expression from [4] 
and create a VHDL-AMS file representing the defective rod. This is then compiled 
through a component level simulator (SMASH 5.2.0™ ||) and can be connected to 
other components. 

The expression for the change in potential energy from the defect can be written 

as 



are the changes to the boundary conditions at the two ends of the beam. 
6. Results 

In order to verify the applicability and accuracy of the defect model (and for that 
matter the ideal Cosserat model) we compare the results of FEA simulations with those 
of the component model. We compare six different slender beams, each a variation of 
a length 150/im, width 6/im and height 15/im beam (see Figure 2). These dimensions 
are quite typical of MEMS devices although we use a shorter beam to ensure that 
the FEA calculations run quickly (about one hour on a 2.80GHz Pentium 4 with 

| SMASH is trademarked to Dolphin Integration 



y(i,o) = / [ n (i,o) . Aw (o,o) m (i,o) . Au (o,o) _ ^,,(0,0)^(1,0)^(0,0) 




- -Au^J^Au^ - -Au^T^Av^ds 
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Case I - Ideal < 




II - As Case I with a blob no 
the fixed end. 

1 SOmierons 



t 

10 in 



Case III - Ideal case where a m 
0.1573ng is attached to the free 
150microns 



Case IV - As Case II with a 
depth 3 at lOOjim. 

150microns 



40 10 10 



Case V - As Case II with 
depth 1.5 at 50fim. 



6microns 

V 

M = 5M B= a m 



Case VI - As Case III with 
depth 1.5. 



6microns 
M=5M „..,,„ 



Figure 2. Diagrams of the Test Cases. Each of these is a variant of Case I and 
the height (15/^m, not pictured) remains unchanged throughout. In all cases we 
consider the lowest energy mode where the motion is in the plane of the diagrams. 




A Case VI - 100, 1.5 

Case V - 50, 1 .5 
Case IV -100, 3 

Case III - mass 

Cass II - bump, no mass 
* Case I - no mass 



(# elements) 1 



Figure 3. Graphs of FEA calculations and their extrapolated values where the 
frequencies are renormalised to a percentage of the Cosserat value. 



1GB RAM for the longest calculations) with a decent accuracy. All simulations are 
performed in FELT [5] and it is worth noting that the solution of the Cosserat model 
takes an imperceptible amount of time on the same computer. 

As a first test we must determine the accuracy of the Cosserat model of 
an ideal beam. For this purpose we must consider the result of a 'perfect' (i.e. 
one with an infinite number of elements in the beam) FEA calculation. For 
this purpose we must extrapolate the results of FEA calculations to infinity. We 
assume that for any FEA calculation, the result will take the form /fea = /oo + 
(5/(numbcr of triangular elements). With this assumption it is possible to calculate 
/oo by extrapolating a number of calculations with varying numbers of elements. This 
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involves making a linear fit to /fea versus l/#elements. Figure 3 demonstrates the 
validity of this method. It is quite apparent from the graphs that the FEA results 
lie close to a straight line for all six cases. Of further assurance is the fact that S 
is positive in all cases. This is as expected since an FEA calculation will necessarily 
be more restrictive than the infinite case and must produce a higher lowest-mode 
frequency. 

We summarise the results of our calculations in Table 1. Considering first the 
two ideal cases (I and III, without and with an attached mass respectively) we see 
that the Cosserat model is a highly accurate model of a beam. In both cases it agrees 
to within 0.5% of the extrapolated FEA result (/<x>)- This is well within the error 
bounds of the extrapolation technique. 



System 


/oo 


/Coss 


% Err 




(GHz) 


(GHz) 




I 


334.2 


335.3 


0.29 


II 


358.3 


358.1 


-0.045 


III 


72.18 


71.85 


-0.46 


IV 


69.16 


70.48 


1.9 


V 


69.28 


69.37 


0.13 


VI 


71.08 


71.22 


0.21 



Table 1. Results of the FEA and Cosserat Simulations. 

Moving on to the defect model we observe that the results of the Cosserat model or 
similarly close to those of the FEA simulations for all cases but IV. The 'poor' accuracy 
of Case IV is easily understood by considering that we have chosen an extreme case 
which goes beyond the expected range of the perturbation technique used. A defect 
this large (half the width of the beam) would not occur in a real MEMS device and 
would almost certainly cause significant problems (such as a full breakage) beyond the 
applicability of any beam model. 

7. Conclusion 

In this paper we have developed a technique for the construction of beam models with 
defects. We apply this technique to the construction of a model of a beam with a 
blob/nick present as well as a model of a beam with random jitter. 

The theoretical framework presented can be effectively applied to new cases for 
the development of new models. This allows us to devise models for defects which 
are presently unobserved (such as kinks in a beam), but which may cause problems in 
future devices, particularly those made of novel materials or with novel applications. 

The developed blob/nick and random jitter models could be used for reliability 
tests of MEMS both at the design and production stage. This is particularly valuable 
in the modelling of high specification MEMS, where an accurate model of faults is 
vital for the design process as costs are high for experiments. 

Tests on the blob/nick model show that the model is a highly accurate tool for 
the simulation of such systems. An accuracy of at worst 0.5% for all but the most 
extreme case (which would not appear in practice) suggests that the model is more 
than satisfactory for simulation of real devices through Component Level Simulations 
or FEA simulations. 



Defective beams in MEMS modelling 



9 



One possible future use of this technique would be to develop a model of a 
curved beam with a defect. At present, curved beams present difficulties for standard 
techniques of ideal beams. The Cosserat model has already been applied to curved 
beam models [11] and application of the aforementioned technique to this work may 
allow accurate modelling of future MEMS devices such as accelerometers. 
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